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Abstract 

Using an innovative damped-Newton method, we report the first calculation 

of many distinct unstable periodic orbits (UPOs) of a large high-dimensional 

extensively chaotic partial differential equation. A majority of the UPOs turn 

out to be spatially localized in that time dependence occurs only on portions of 

the spatial domain. With a particular weighting of 127 UPOs, the Lyapunov 

fractal dimension D = 8.8 can be estimated with a relative error of 2%. We 

discuss the implications of these spatially localized UPOs for understanding 

and controlling spatiotemporal chaos. 
47.27.Cn, 05.45.+b, 05.70.Ln, 82.40.Bj 
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Much recent research on sustained nonequihbrium systems has concerned using unstable 
periodic orbits (UPOs) to characterize and to control chaos The starting point for 

this research has been the observation that a dense set of unstable periodic orbits (UPOs) 
is associated with a strange attractor. 0. For low-dimensional chaos with a symbolic 
dynamics (a unique labeling of the UPOs in the dense set |]^), researchers have shown that 
knowledge of a small number of UPOs can be used to improve forecasting 0] and to estimate 
dynamical invariants of the strange attractor such as its fractal dimensions and Lyapunov 
exponents In some low-dimensional chaotic experimental systems, enough UPOs have 
been determined from empirical data to apply this formalism successfully 0. An exciting 
related advance has been the discovery that these low-dimensional chaotic systems can be 
controlled by stabilizing particular members of the dense set of UPOs with weak time- 
dependent perturbations of a system parameter P]. 

An important question is whether the above results generalize and remain useful for 
high-dimensional dynamical systems, especially for large nonequihbrium systems that are 
extensively chaotic 0,[T^. We are particularly interested in two potential applications. One 
is to characterize extensive chaos in terms of some finite number of UPOs 0]. Large nonequi- 
hbrium chaotic systems are difficult to analyze because of long transient times (which, in 
some cases, increase exponentially with system size and also because it is often not 
known in advance how long a system should be observed in a non-transient regime to obtain 
good statistics. An attractive feature of UPOs associated with a strange attractor is that 
they can be computed within a finite time interval and so their analysis replaces difficult 
questions of transients and observation times with questions of mathematical and numerical 
convergence, e.g., the number of UPOs needed to approximate the invariant measure of an 
attractor to some desired accuracy and whether sufficiently many UPOs can be calculated 
numerically. A second potential application is to improve control of extensively chaotic sys- 
tems. If the magnitude of a control perturbation is to remain small as a chaotic system 
becomes large, the system parameter will need to be varied at spatially distributed control 
points |]12|. An improved understanding of the spatial structure of UPOs, of the distribu- 



tion of their periods T, and of their stabihty should aid the development of high-dimensional 
spatiotemporal control algorithms by suggesting the number and location of control points 
for a particular UPO and for a particular system parameter. 

In this Letter, we take a significant step towards understanding the relation of the dense 
set of UPOs to high-dimensional spatiotemporal chaos by reporting the first calculation of 
many (over 100) distinct UPOs for a high-dimensional {D = 8.8) driven-dissipative system 
partial differential equation (pde) in an extensively chaotic regime [T^. This calculation 
represents two achievements. One is numerical, that a simple modification of a Newton 



algorithm by the addition of damping |jTj] greatly increases the likelihood of convergence 
and so makes practical the computation of many UPOs from the nonlinear equations that 
they satisfy. The second achievement is several discoveries in nonequilibrium physics made 
possible by this numerical method: that over 100 distinct UPOs can be calculated in a large 
chaotic pde and therefore their properties can be studied for the first time; that most of 
these UPOs turn out to be spatially localized (as discussed below); that about 100 UPOs are 
already sufficient to estimate the fractal dimension of the high-dimensional attractor to two 
significant digits; and that certain mean quantities such as a time-averaged spatial pattern 
can be estimated from a knowledge of the UPOs. 

Especially interesting is the result that most of the UPOs have dynamics that are spatially 
localized in that the variance of temporal fluctuations is substantial only on portions of the 
spatial domain (see Figs. |l] and below). We speculate that this will be a general feature 
of UPOs for extensively chaotic systems. This localization, in fact, helps to explain how 
the spatiotemporal disorder of an extensively chaotic system (which can be interpreted as 
many small dynamically-independent subsystems [|1^]) can be consistent with the global 



spatiotemporal coherence of each UPO in the dense set. The temporal disorder can be 
understood as an irregular wandering of a chaotic orbit near many UPOs (the case for low- 
dimensional chaos) while the spatial disorder can now be better understood as a complicated 
weighting of localized spatially-irregular dynamics associated with the UPOs. The spatially 
localized dynamics also has implications for distributed control algorithms since control 



points placed in the weakly time-independent spatial regions of a particular UPO will likely 
not aid in stabilizing that UPO. It is then important to know the spatial structure of various 
UPOs before attempting to control them. 

Our calculations were carried out for one of the simplest models of extensive chaos, the 
one-dimensional Kuramoto-Sivashinsky (KS) equation ||^ 

dtu = —u dxU — d^u — d^u, x G [0, L], (1) 

for a field u{t,x) that lives on an interval of length L and that satisfies "rigid" boundary 
conditions u = dxU = at x = and x = L. For system sizes L > 50, Manneville has 
shown that typical initial conditions evolve towards a chaotic attractor that is extensive in 



that the Lyapunov fractal dimension D increases linearly with L ||T5|. In our calculations, 
we chose a fixed length L = 50 and spatial resolution Ax = 0.5 for which the Lyapunov 
fractal dimension was D = 8.8 and there were 4 positive, 1 zero, and 94 negative Lyapunov 
exponents . The system size L = 50 was just large enough to be in the extensively chaotic 
regime and yet small enough that the numerical calculations were manageable with available 
resources and algorithms. The reason for emphasizing the extensively chaotic regime was 
that certain details of the spatiotemporal chaos and of the UPOs could then be expected to 
be insensitive to the value of the system length L. 



UPOs of Eq. (|T]) were calculated numerically by using a damped-Newton method |T4 



with simple shooting [|T^, pages 120-122] to solve a set of nonlinear equations F(X) = for 
unknowns X = (T, Uq) which we define as follows. Given a uniformly spaced mesh Xj = iAx 
on the interval [0, L] with spatial resolution Ax = L/{N + 1), we denote the field value at 
the ith mesh point at time t by Ui{t) = u{xi, t) for i = 0, ■ ■ ■ , -|- 1. Then the (A^ -|- 1)- 
dimensional vector of unknowns X consists of the UPO's period T and of a A^- dimensional 
point Uq = (ui(0),- ■■,U]\f{0)) that lies on the UPO. The values Uq = u{t,0) and mat+i = 
u{t, L) are zero by the boundary conditions and do not need to be solved for. 

The corresponding (A^ -|- l)-dimensional nonlinear equations F(X) = for X are de- 
termined by trying to find a closed orbit \Jt of duration T starting from Uq, i.e., such 



that Ut — Uo = 0. These nonhnear equations are derived and solved as follows. First, 
Eq. (|1]) is reduced to a iV-dimensional set of first-order autonomous odes, 



dV/dt=G{V), U(0) = Uo. 



(2) 



by replacing the spatial derivatives in Eq. (|T|) and in the boundary condition d^u = with 
second-order accurate finite-difference approximations. The vector Ut at time T is then 
found by numerical integration of Eq. (|^). Second, a N x N matrix at time T is 
obtained by integrating the N"^ linear variational equations [O 



dM/dt = JM, M(0) = I, 



(3) 



together with Eq. (^, where J(t) = dG/d\J denotes the N x N Jacobian matrix and 
I denotes the N x N identity matrix. (The matrix M(t) is the propagator that evolves 
infinitesimal perturbations from time t = to time t.) Given the vector Ut and matrix My, 
a Newton correction 5X = {6T, 6\J) for updating the vector X is found by solving the 
following 2x2 block matrix equation [p!^ : 



G(Uo, 
G(Ut) Mt-I 





Uo-U 



T 



(4) 



where the symbol f denotes the matrix transpose. A single Newton step is then performed 
by updating the present values for the unknowns, X X + (5X, and such steps are repeated 
until the residuals and corrections are sufficiently small 



|Ur-Uo||oo < lO-^U, 



II OO) 



and 



|(5X|L < lO^^llX, 



oil OO) 



(5) 



in the infinity norm ||X||oo = maxj |Xj|. Each Newton step requires the time-integration of 
a set of + A^^ odes over time T plus the solution of the linear equations Eq. (^). For the 
calculations reported below, we used time-splitting methods for Eqs. (0) and (|^) that were 
first-order accurate in time and second-order accurate in space with a constant time step 
At = 0.005. Since the block matrix — I in Eq. (^) is typically dense, PL[/-factorization 
was used to solve for the correction; the time to carry out the linear algebra was much larger 



the time to integrate Eqs. and @ over a period T. Most UPOs that we calculated were 
robust to modest changes in spatiotemporal resolution and the Newton convergence criteria. 

For high-dimensional Newton methods, it is essential to have a good starting guess Xq 
since Newton methods are guaranteed to converge only locally. We initially tried to find 
a good initial guess (T, Uq) by searching for approximate recurrences [^,0 of chaotic time 
series = U(iAt) in the high-dimensional numerical phase space of Eq. (|^). This turned 
out to be impractical, e.g., for L = 50, integration times of at least 10*^ time units were 
needed to find a single approximate recurrence of period T = 11.6 within a ball of rather 
large radius 0.1, ||Ut — Uo||oo < 0.1 [|l^. Further, these approximate recurrences did not 
turn out to be close to any of the UPOs that we calculated using the more rigorous Newton 



method [|T9|. Since these results suggest that no approximate recurrence is close to a UPO 
of Eq. (|l]), we then tried to choose an initial state by choosing a positive random number T 
for the period and an initial vector Uq from a point on the chaotic attractor. Using this 
procedure, the convergence criteria Eq. always failed within the specified maximum 
number of 200 Newton iterations. 



Convergence was finally obtained with a damped-Newton method ||T3, in which only a 



fraction a < 1 of a Newton correction (5X was added to update the unknowns, X X+a5X. 



We used a particular damping method known as the Armijo rule flj] which expresses the 
damping factor in the form a = 2~K During each Newton iteration, successive values of 
the integer j = 0, 1, ■ ■ ■ were tested until a value j was found such that the new residual 
||F(X + 2^-'(5X)|| was smaller by a factor 1 — q;/2 than the previous residual ||F(X)|| in 
the Euclidean norm. With the Armijo rule, Newton's method converged with the criteria 
Eq. for roughly five percent of all initial guesses (T, Uq) consisting of a positive random 
number T and some point Uq on the chaotic attractor. If the damping parameter a fell 
below the value 2~^, it was more efficient to terminate the Newton iteration and start over 
with a new guess (T, Uq). The specific value of j = 6 was chosen as a trade off between 
convergence and computer time, as tiny fractions of the correction can lead to an excessive 
amount of time spent on a single guess. 



We now discuss the properties of the UPOs calculated with the above numerical methods. 
Using the Armijo rule and the Newton method discussed above, 127 distinct UPOs were 



found from 5000 guess Xq = (To,Uo) [T^. UPOs with periods shorter than 8 could not be 
found while the Armijo-Newton algorithm failed to converge for UPOs with periods larger 
than 42, probably because the simple shooting method becomes unstable for periods larger 
than about this value. 

As shown qualitatively in Fig. |l] and more quantitatively in Fig. ^, most of the UPOs are 
localized in that their time variation is substantial only in isolated portions of the domain. 
Fig. |l](A) shows a state for which the temporal variation is localized in two small spatial 
intervals, one of which is adjacent to the right boundary. Fig. [l|(B) shows an interior state 
whose dynamics is localized to a single interval well away from the boundaries. Fig. |l](C) 
shows an extended UPO in which the temporal variation occurs throughout the spatial 



domain. The corresponding mean and variance patterns [|19| are given in Fig. 0(a)-(f). In 



Fig. 0(a), (c), and (e), the mean pattern is nonzero throughout the domain which holds for 
the other 124 UPOs as well. Fig. 0(b) and (d) indicate more clearly the localization of the 
dynamics, which is evidently uncorrelated with the mean pattern. The variance decreases 
three to five orders of magnitude outside the regions of substantial variation. 

From an explicit knowledge of the space-time evolution of the 127 distinct UPOs, the 
time-averaged mean and the variance of the spatiotemporal chaotic solution of KS equation 
Eq. (|l]) can be estimated. A simple averaging < m{x) >upo of 127 mean patterns m{x) 
for the UPOs yields the solid curve in Fig. 0(g), which should be compared with the dashed 
curve obtained by averaging an extensively chaotic field u{t,x) over 10^ time units. The 
relative error in the infinity norm between the patterns is 24% and is substantially better 
near the boundaries. Strikingly, < m{x) >upo has the same qualitative structure as the 
mean pattern of the attractor < u{x) > using only a moderate number of low-period UPOs. 
In contrast, an average of the 127 variance patterns (solid line in Fig. 0(h)) does not agree 
as well (a relative error of 46%) with the variance of the extensively chaotic field u{t,x) 
averaged over 10^ time units but still reflects some of the qualitative features. These results 



suggest that a knowledge of UPOs may help to understand the interesting mean patterns 
and localized dynamics observed recently in Faraday wave experiments [ pO[| . 

Fig. 1^ shows how the extent of localization and instability depends on the period T. To 
characterize the degree of localization of the UPOs, a localization number was defined as the 
fraction of the interval [0, L] for which the variance v{x) was smaller than 0.05 (the results 
were not sensitive to the choice of this cutoff). Fig. ^(a) shows the localization number 
versus the period T for all 127 UPOs. Although there is scatter in the points, we can draw 
a few conclusions: that UPOs with period less than about 8 don't exist; that shorter period 
UPOs are more strongly localized; and that there can be many UPOs of approximately the 
same period (say T = 14) and these can vary substantially in their localization. The outlier 
point with a localization of 0.6 for period T = 30 suggests that the right side of this curve 
is incomplete, that there are localized UPOs of higher period but our numerical algorithms 
were not able to find them. In Fig. 0(b), we summarize the instability of all 127 UPOs as 
a function of their period T. The instability of a UPO is defined as the sum X]+ of all 
positive transverse Lyapunov exponents which are defined by A = log(|m|)/T, where m is a 
Floquet multipher, i.e., an eigenvalue of the propagator matrix in Eq. (0) evaluated for 
a UPO of period T |T^. Fig. ^(b) shows that, on average, smaller period UPOs are more 



unstable. Again the right side of this curve is incomplete and it would be interesting to 
extend the data using more sophisticated numerical algorithms. 

Using the data in Fig. ^(b) one can estimate the fractal dimension D of the chaotic 
attractor. First, a fractal dimension is formerly associated with each UPO by generalizing 



the Kaplan- Yorke formula [jT^ to UPOs in terms of their transverse Lyapunov exponents; we 
find dimensions ranging from 6 to 12 for the 127 UPOs. The fractal dimension of the chaotic 
attractor can be viewed as an average over the fractal dimensions of the UPOs which are 
associated with the attractor. To obtain the fractal dimension over attractor, each UPO's 
dimension is summed with a weighting 1/(I]+A) (less unstable UPOs are weighted more 
heavily) which gives an estimate of -D = 9.0 to the fractal dimension of the attractor. This 
is in good agreement with the Lyapunov fractal dimension D = 8.8 ±0.1 calculated directly 
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from the Lyapunov exponents of the spatiotemporal chaotic solution of Eq. (|T|) |T^. The 



convergence to the Lyapunov dimension of the estimate based on UPOs is statistical in that 
the error decreases approximately as 1/\/N where is the number of UPOs contributing 
to the weighted sum. Other weightings of the UPOs were tried [0,^ but found not to give 
results as satisfactory as the 1/(1]'^+) weighting. 

In conclusion, we have used an innovative numerical method to calculate for the first time 
many UPOs associated with a large high-dimensional pde, Eq. ([l|). An important numerical 
insight was the use of damping to increase the likelihood of convergence of an otherwise 
straightforward Newton method. Near-recurrences in the time-series data as long as 10^ 
time units were found not to correspond to any computed UPOs. From the 127 distinct 
UPOs found using the damped-Newton method, we could estimate the fractal dimension of 
the attractor to be 9.0 ±0.1 compared to the actual value of 8.8 ±0.1. These UPOs were also 
used to predict successfully the qualitative features of the time-averaged mean pattern and 
the variance of the chaotic attractor. A new insight is that the UPOs are typically localized 
in space. This localization suggests a new way to think about the dynamically independent 
subsystems within extensive chaos. The localization also has important implications for the 
control of large chaotic systems using distributed sets of control points. 

The present calculations raise two interesting questions for future work. One is to improve 



the numerical algorithms, primarily by using multishooting methods [^, so that UPOs of 
partial differential equations in two- and three-space dimensions can be calculated. It will 
then be quite interesting to determine whether there are systematic patterns of defects asso- 
ciated with UPOs as their period increases from their smallest value. The second question 
is to understand the implications of localized UPOs for spatially distributed control, e.g., 
by the OGY or delayed-feedback methods. 

We thank L. Bunimovich, S. Newhouse and M. Strain for useful discussions. This 
work was supported by a DOE Computational Science Graduate Fellowship, by NSF 
grants NSF-DMS-93-07893 and NSF-CDA-92123483-04, and by DOE grant DOE-DE-FG05- 
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FIGURES 

FIG. 1. Density plots of three representative UPOs u{t,x) calculated by applying a 
damped-Newton method to the Kuramoto-Sivashinsky equation Eq. (||) in a spatial domain of 
length L = 50. The horizontal axis is space and the vertical axis spans a time interval of 35 time 
units, (a) A UFO of period T = 9.9 with dynamics localized near the right boundary; (b) A UFO 
of period T = 10.7 with dynamics localized in interior of the interval; (c) An extended UFO of 
period T = 23.4. The greyscales represent amplitude variations between about 3 and -3. 

FIG. 2. Time-averaged mean patterns m{x) =< u{x,t) > and variance pat- 
terns v{x) =< — m(x))^ > for the three representative UFOs of Fig. |l]. (a) and (b): 
for the UFO with dynamics localized near a boundary, (c) and (d): for the UFO with dynamics 
localized away from boundaries, (e) and (f): for the extended UFO. (g) and (h): mean and 
variance patterns (solid lines) averaged over all 127 distinct UFOs. For comparison, the dashed 
lines give the corresponding mean and variance patterns obtained from an integration of a chaotic 
solution over 10^ time units. 

FIG. 3. (a) Localization (fraction of the spatial domain that has variance v(x) below 0.05) 
versus period T of 127 distinct UFOs calculated for the KS equation (Eq. (||)) in an extensively 
chaotic regime with system size L = 50. (b) Degree of instability as measured by the sum ^_|_ A 
of positive transverse Lyapunov exponents versus the period T for all 127 UFOs. 
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